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ABSTRACT 

A commonly cited drawback of SPH is the introduction of spurious shear viscosity 
by the artificial viscosity term in situations involving rotation. Existing approaches 
for quantifying its effect include approximate analytic formulae and disc-averaged be- 
haviour in specific ring-spreading simulations, based on the kinematic effects produced 
by the artificial viscosity. These methods have disadvantages, in that they typically 
are applicable to a very small range of physical scenarios, have a large number of 
simplifying assumptions, and often are tied to specific SPH formulations which do not 
include corrective (e.g., Balsara) or time-dependent artificial viscosity terms. In this 
study we have developed a simple, generally applicable and practical technique for 
evaluating the local effect of artificial viscosity directly from the creation of specific 
entropy for each SPH particle. This local approach is simple and quick to implement, 
and it allows a detailed characterization of viscous effects as a function of position. 
Several advantages of this method are discussed, including its ease in evaluation, its 
greater accuracy and its broad applicability. In order to compare this new method 
with existing ones, simple disc flow examples are used. Even in these basic cases, the 
very roughly approximate nature of the previous methods is shown. Our local method 
provides a detailed description of the effects of the artificial viscosity throughout the 
disc, even for extended examples which implement Balsara corrections. As a further 
use of this approach, explicit dependencies of the effective viscosity in terms of SPH 
and flow parameters are estimated from the example cases. In an appendix, a method 
for the initial placement of SPH particles is discussed which is very effective in reducing 
numerical fluctuations. 

Key vifords: Hydrodynamic models - Viscosity - Accretion and accretion discs 



1 INTRODUCTION 

Continued advances in computing power and availability 
have led to an increasingly important role for numerical 
simulations in understanding hydrodynamical phenomena, 
and astrophysics provides a particular range of interesting 
applications. High-resolution, multidimensional studies fre- 
quently include an ambitious amount of physics in order to 
reproduce the dynamics of a system with sufficient accuracy. 
However, numerical artefacts and limitations are necessarily 
present and must be understood as fully as possible to ap- 
preciate the actual physics in the systems being represented. 
Non-physical features enter ab initio into all hydro- 
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dynamic simulations simply because of the discretisation 
of a 'continuous' fluid (which is itself discrete, of course, 
but on very small scales). The two types of method de- 
veloped for this task have very different structures: La- 
grangian methods follow moving fluid elements, while Eule- 
rian methods compute flow properties on a grid of points; ad- 
ditionally, there are 'hybrid' approaches, the so-called ALE 
(arbitrary L agrangian-Eulerian) schem es, e.g. the BETHE- 
hydro code (|Murphv fc Burrows! 120081 '). In this paper, we 
focus on the Lagrangian smoothed particle hydrodynamics 
(SPH) method, which ha s been widely used for modelling 
astro physical phenomena (|Gingold fc Monaghan|[T977l : ILucvI 
Il977h . SPH has both advantages and disadvantages with re- 
spect to other numerical techniques. One main advantage is 
that it automatically adapts to following dynamic flows and 
arbitrary geometries, without the need for mesh reflnement 
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or other readjustment techniques required with Eulerian- 
based codes. Furthermore, the energy equation is solved in 
the local comoving frame of each fluid element, giving a di- 
rect implementation of the first law of thermodynamics. 

In modelling fluids by means of SPH particles, most 
SPH codes have included an artificial viscosity term, both for 
making the ensemble of discrete particles behave more like 
a continuum (in continuous regions of the flow) and for han- 
dling shock discontinuities which may arise during the course 
of a simulation. Methods for shock-handling with Ricmann 
solvers have also been im plemented in SPH (|Inutsuka,2002i : 
ICha fc Whitworthll2003l l: these approaches have both ad- 
vantages and disadvantages. Here, we focus on the artiflcial 
viscosity methods as used in standard SPH formulations. As 
well as having beneflcial effects, artificial viscosity can also 
introduce unwanted numerical artefacts, among which spu- 
rious shearing torques in r otating flows are par ticularly trou- 
blesome (see, for example. iFlebbe at al.l l|l994h ). Techniques 
such as the Balsara correction (described in the next section) 
have been developed for counteracting this. Typically, SPH 
prescriptions for artificial viscosity are a composite of both 
bulk and shear viscosity components, but act as a source 
of dissipation primarily related to shear when modelling ro- 
tating fluid flows. It is important to have some quantitative 
measure of how large the effects of this are, whether one 
is dealing with problems having no intrinsic shear viscosity 
or whether there is also a real, physical shear viscosity; in 
either case, any numerical viscosity should be suitably neg- 
ligible (unless one is trying to use an arti ficial viscosity to 
actually model a physical viscosity, see e.g.. lLodato &: Pried 
l|2010l )). We focus here on two standard artificial viscosity 
prescriptions, with and without the Balsara correction. 

The main objective of this paper is to advocate moni- 
toring local entropy generation as a means for quantifying 
the residual effects of artificial viscosity in regular flows not 
involving shocks. Our method is completely general, but 
we test it out by applying it to frequently-used idealized 
test problems involving isothermal discs and compare its 
performance on these with that of other alternative meth- 
ods. Two main types of approach have been used previ- 
ously, both of which measure the kinematic effect of the 
artiflcial viscosity and relate it to an effective viscosity co- 
efficie nt, v. First ly, there are the empirical 'ring-spreading' 
tests (lLustlll9 52': 'Lvnden-B eU fc Pringlelll974l : iFlebbe et al] 
[l99i;|Murr3l996; Spcith fc Riffertlll999[ l. These are based 
on the analytic relations for the kinematic contribution of 
the artificial viscosity to disc-averaged behaviour, with as- 
sumptions of an isothermal equation of state and a small 
viscosity (so as not to generate large radial flow) that is 
constant across the disc. In practice, these tests can only be 
performed in two dimensions, as the thermal pressure of the 
fluid is zeroed to isolate the viscosity term in the equations 
of motion. These methods also require post-simulation flt- 
ting of averaged results. Secondly, analytic approximations 
have been develop e d relating v to SPH and disc parame- 
ters llMurravllTggg : iLodato fc Pried l201Ct ). These relations 
are generally derived for specific formulations of SPH arti- 
ficial viscosity with the quadratic term set to zero. Various 
additional assumptions include having constant smoothing 
lengths (in the continuum limit) and particular SPH par- 
ticle kernel shapes, as well as considerations that the vis- 
cosity is active for both receding and approaching particles. 



Both of these kinds of existing approach involve restrictive 
assumptions and neither may be applied to artificial viscos- 
ity prescriptions that include often-used 'corrective' terms, 
such as the Balsara method or time-dependent approaches 
(discussed below), which aim to reduce 'excess' shearing. 

The technique presented here involves evaluating the 
effect of artificial viscosity in rotating flows directly from lo- 
cal entropy production (i.e. using the energy equation rather 
than the equation of motion). It does not assume any par- 
ticular equation of state or rotation proflle and also does 
not assume constancy of v or of smoothing lengths, or any 
particular artificial viscosity prescription. Moreover, it re- 
quires no special simulation setup and is directly calculable 
from standard SPH particle quantities. Impetus for this ap- 
proach arises from the nature of the SPH method itself and 
its direct application of the first law of thermodynamics, 
noting that a clear understanding of entropy production is 
of great importance in numerical simulations. The technique 
is simple and quick to implement, and it allows a detailed 
characterization of the effects of the artificial viscosity as a 
function of position. Several advantages of this method over 
existing ones are discussed, including its ease in evaluation, 
its direct interpretation, its greater accuracy and its broad 
applicability, e.g., to arbitary rotation profiles and equations 
of state. 

For all of the simulations described here, we have 
utilised the pu blic-release version of the Gadget-2 code 
(|Springell I2OO5V While several artificial viscosity formula- 
tions of varying sophistication have been developed over the 
years, in this paper we focus mainly on the method of es- 
timating the effects of this viscosity and therefore utilize 
only the relatively simple form that is standard in Gadget-2 
(with and without the Balsara treatment, described below). 
However, the method itself is independent of viscosity pre- 
scription, and our aim is to apply it to a comparison of 
several other formulations in a forthcoming paper. Here, we 
begin by briefly reviewing in ^the form of SPH equations 
used and considerations concerning the artiflcial viscosity; 
in §31 the ring-spreading test is described; in §3] both the 
motivation for and the general formulation of our local cri- 
terion for measuring the effective viscosity in rotating flows 
are described; in ^ example applications of our method 
in thin discs are given, and various simulations performed, 
for which the results of the local measures of viscosity are 
compared with those of the ring-spreading tests and existing 
analytic formulations; these results are then discussed in ^ 
Additionally, in Appendix |^ we discuss the SPH particle 
setup used for our initial conditions, which is presented as 
a simple improvement to aid in reducing numerical effects 
and increasing the efficiency of convergence. 



2 SPH, USING GADGET-2 

In SPH a continuous fiuid is sampled at a finite number 
of points, and discretised versions of the standard hydrody- 
namic equations of continuity, momentum (Euler) and the 
first law of thermodynamics (energy equation) are evolved, 
given respectively in their simplest form as: 



= 5^P + PV.V 



(1) 
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= -iv + ivP 

at p 



= 



P d 



dt 



(2) 
(3) 



written with Lagrangian derivatives, d/dt = d/dt + v ■ 
V, and with mass density, p; velocity, v; thermal pres- 
sure, P; and specific internal energy, e. The exact form of 
SPH implementations has continued to be refined over the 
decades; in this work, we have used the publicly- available 
SPH code, Gadget-2, the formulation of which is derived 
from a fluid Lagrangian and, importantly, directly con- 
serves linear and a ng ular m omentum, entropy and energy 
jSpringel fc Hernauistll2002h . The Gadget-2 SPH equations 
are presented here, with emphasis on the necessity for in- 
cluding extra, purely numerical terms in the implementa- 
tion of the continuum equations above. For further details, 
we refer the rea der to mu ch more com p rehen sive discussions 
elsewhere, e.g. (|MonaghanH2005l : |Priceli2010l ) . 

The fluid is taken to have a polytropic equation of state 
(EOS), 



P = K{a)p-' = (^-l)pe, 



(4) 



wher e s is the specific entrop y. Following standard usage 
fe.g.. lLandau fc Lifshitj (|l959l ')'). a polytropic gas is defined 
as having constant specific heats, and therefore the related 
adiabatic index, 7, is constant. While some authors restrict 
K to being constant as well (and refer to it as the poly- 
tropic constant), here it is allowed to vary, and in the fol- 
lowing K is referred to as the entropic function. In Gadget-2, 
this K is evolved rather than e, giving the strict conserva- 
tion properties of both e nergy and entropy mentioned above 
l|Springel &: Hernauistl2002 ). In the following, the subscripts 
{a, b} will be used only to distinguish individual fiuid ele- 
ments, represented in SPH as finite volume pseudo-particles, 
and not to signify vector components. 

First, the mass density at the location of a given parti- 
cle labelled as 'a' and with spatial coordinate vector, r^, is 
calculated at each timestep by direct summation using 



pa = ^ nitWa . 



(5) 



The summation is performed over A'n 'nearest ne ighbour' 
parti cles labelled as (typically ~ 50 (Navarro fc White! 
I1993I )). each with mass, nib, and contained within a ra- 
dius given by the smoothing length, ha- The kernel, Wa = 
VK(|rai,|, /la), is strongly peaked and differentiable (a cubic 
spline, in Gadget-2) and with |rai,| = |ra — ri,| < ha- The 
smoothing length is calculated at each timestep such that 
the kernel volume of a particle contains a constant mass (for 
numerical stability), according to the following relation: 



— n-apa = N-n_m — coust. 



(6) 



where m is the average particle mass. Within the fluid, the 
resolution is essentially determined by the magnitude of ha, 
over which discontinuities are smoothed. 

To evolve the particle velocity, the non-viscous Euler 
equation, including a gravitational potential, (j), is repre- 
sented as: 



dva' 
. dt . 



b 



Pb 



PS 



where 



fa 



1 + 



hg dpa 
3pa dha 



(8) 



The factors, fa, arise from the Lagrangian derivation and the 
set of constraints on coordinates provided by Eq. [B] they 
account dir ectly for the variation in s moothing lengths in 
the system jSpringel fc Hernauist|[2002l ). Self-gravity among 
particles is calculated efficiently using a tree algorithm (a hi- 
erarchical multipole expansion) (|Barnes fc Hut|[l986l ). Addi- 
tionally, in some cases we include the acceleration due to a 
(Newtonian) central object located at the origin of the sim- 
ulation's coordinates, V0a ~ —GAlcia/ra- 

Gadget-2 uses an artificial viscosity of the form 



Haft = 



^viCabWab — Swlb)/ 2 Pab - - 

(9) 





for Vai, ■ Tab < 0, 
for -Vab ■ Tab ^ 0, 



with Vab = Va - Vb and Wab = {'Vab ■ rab)/irab| being the 
relative velocity between particles in vector and scalar form, 
respectively; pab = (Pa + pb)/^, the average density; Ca — 
i'yPa/ pa)^^^ , the local soundspeed, with average Cab = {ca + 
Cb)/2; and Ov ~ 1, a free parameter for the strength of the 
viscosity. In practice, H^b acts similarly to a pressure term 
in the Euler equation, and the related acceleration is: 



dva 



-^] = -I]"^6n,bVVKab, 



(10) 



where Wab = {Wa + Wb)/2. 

The Hab formulation possesses two terms to mimic 
naturally-occuring, dissipative processes. The first (linear in 
Wab) provides both the bulk and shear viscosity of the con- 
verging particles, and the second (quadratic in Wab) func- 
tions as a von Neumann-Richtmyer artificial viscosity for 
shock handling and for spreading shock discontinuities over 
the smoothing length of supersonic particles. The quadra- 
trie term also prevents interparticle penetration. It should 
be noted that this single artificial viscosity provides both 
shear and bulk viscosities, which would be given by sepa- 
rate terms in the physical, Navier-Stokes description (direct 
implementations of the latter are discussed briefiy below). In 
many SPH prescriptions the linear and quadratic terms are 
scaled by separate free parameters, Ov and /3v, respectively, 
each with a range of 'typical' values, but we here follow the 
practice of setting /3 = 3a (unless otherwise stated). An 
in-depth discussion of the relati ve scaling of the linear and 
quadratic terms is provided by (|Monaghanl 1 19971 ) . 

The entropic function, K (defined in Eq.[3|, changes, in 
general, in the presence of viscosity: 



dKa 
dt 



7 ■ 



pa 



bVab 



VM^a 



(11) 



Cooling in optically thin material may be included in a 
straightforward manner by subtracting the appropriate en- 
tropy loss at the end of a timestep. The inclusion of accurate 
radiative transport in optically thick material is a separate 
(and highly non-trivial) task, which we do not consider here. 

Finally, to limit spurious angular momentum transfer 
which arises in shear fiow due to the form of the artificial 
viscosity, the simple Balsara co rrection is commonly utilised 
l|Balsaralll995l : ISteinmetj|l996l '). The aim of this factor is to 
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remove the effect of the artificial viscosity when there is a 
pure shearing motion and to have it acting only when there 
is a compression. An estimate of the relative amount of local 
shear is made from the curl and divergence of the particle 
velocities: 

|V XVU + |V- VU + (0.0001 Ca//la) ' ^ ' 

where the third term in the denominator prevents the quan- 
tity from diverging. The average value of this quantity for 
two interacting particles is used as a multiplicative fac- 
tor in front of Ilab, producing the modified form Il^f, = 
flat {ga + gb)/2, which is then used instead of Hat in Egs.fTUl 
and llll In the following, we examine the effect of the artifi- 
cial viscosity both with and without the Balsara correction 
(the former case being the default). 

Egs.lSirTlandfTTl together with the EOS and the smooth- 
ing length condition in Eq. (6] form a complete set of equa- 
tions for evolving the physical, hydrodynamic quantities in 
numerical simulations. However, as part of the process of 
discretisation, a number of additional variables and param- 
eters have been included {h, W, etc.). While these are used 
to model realistic properties of the discretised system, such 
as continuity and shocks, they also introduce purely numer- 
ical features by controlling resolution, the 'spread' of shocks 
and the stability of interpolations and summations. In the 
past decades, much work has been done to determine reli- 
able values and forms for these parameters and to reduce 
associated numerical artefacts. Generally, a number of stan- 
dard problems, in various dimensions and geometries, are 
used to test the behaviour of a given code and the associ- 
ated parameters. 

Another purely numerical consideration in using SPH 
is the initia l placement of the partic l es representing p oints 
in the flow (jlmaeda fc Inutsukal[200i : lMonaghanll200d ;). For 
instance, regularities in particle configurations may lead 
to the propagation of artificial structures in the simula- 
tion along preferred directions, and even randomised place- 
ments may produce numerical features, e.g., due to local 
overdensities. For each simulation presented here, the ini- 
tial particle setup was created using an algorithm which 
maps a non-regular but constant number density profile 
of (equal-mass) particles into the arbitrary number-density 
(and trivially, mass-density) profile of the desired model. 
This method has also been applied in other simulations 
l|Tavlor. Miller fc Podsiadlowskil I2OIII ) and is described in 
Appendix El 



3 EXISTING (KINEMATIC) 
APPROXIMATIONS OF f 

The standard ring-spreading test starts from an initial 5- 
function ring of matter in circular Kcplcrian motio n around 
a gra vitating point mass yLust 1952; Lyndcn-BcU fc Pringld 
[l97J). If shear viscosity is present, this ring will proceed to 
spread out into a disc, with the rate of spreading dependent 
on the magnitude of the viscosity. The matter is taken here 
to behave isothermally, with the temperature and viscosity 
being sufficiently small so that the spreading is slow and 
the thin disc is essentially Keplerian (with its height, H, 
being small compared with the radial coordinate, r, at any 



point) (|Pringle|[T98ll '). With certain additional simplifying 
assumptions (negligible self-gravity for the disc material and 
constant viscosity coefficient, v), an analytic solution can be 
obtained for the ring-spread ing (|Flebbe et al.|[l993 : lMurravl 
1 19961 : ISpeith fc Rifer^ll999l ') . 

The exact analytic results provide a means for com- 
parison for an SPH simulation of an equivalent configu- 
ration with no viscosity apart from that of the numerical 
scheme, i^sph. From the analytic results, an effective value 
for the shear component of the artificial viscosity can then 
be estimated, if one assumes that vsph is roughly constant 
throughout the disc (a somewhat doubtful assumption) and 
that the corresponding bulk component has negligible ef- 
fect (a very good assumption under the circumstances en- 
visaged). The value obtained can then be translated into 
an effective Shakura-Sunyaev a coefficient [25], which rep- 
resents the viscosity according to 

u — acs H, (13) 

where H is the disc height, and Cs is the sound speed. 

The analytic solution gives an expression for the sur- 
face density of the ring, Er, as a function of position and 
time. For the conditions mentioned above, the radial veloc- 
ity Vr (positive for 'outer' material and negative for 'inner' 
material) is always small, with Vr ^ v jr <C by scaling ar- 
guments (|Frank. King fc Rainel l2"002h . The surface density 
is then given by 

E.(.,.)^^exp(^)/,,,(^^), (14) 

where r is the radius in the equatorial plane (with r = 1 
initially); r is a time parameter (r = 121/1); v is the kine- 
matic shear viscosity coefficient; A/d is the disc mass (taken 
as the mass unit, i.e. Md = 1), and /1/4 is a modified Bessel 
func tion. A detailed deriv a tion o f this formula is given by, 
e.g. , iFrank. King fc Raind (|2002l ) . 

In our work reported below, the corresponding numer- 
ical simulations begin a short time after the axisymmetric 
ring has spread to a finite thickness (at ro = 0.01) using 
the analytic formula Eq. [TJ] to provide the initial surface 
density profile. This same formula is also used for match- 
ing against the simulation properties at subsequent times 
(doing so has been justified in previous numerical studies 
by the general proximity of the surface density evolution 
to this expression). Self-gravity of the SPH particles is in- 
cluded in simulations, although it is assumed that the ac- 
celeration due to the central point mass dominates (we take 
Mc/A'/d = 1000). SPH particles are removed from the sys- 
tem at an inner boundary with small radius, r = 0.1. 

As mentioned above, in addition to the ring-spreading 
test there are also some approximate analytic formulations 
for estimating the effective shear viscosity of the disc from 
kinematic considerations: that of iMurravi (|199S ) is 

vm = avCs/i/16 , (15) 

while that of iLodato fc Pried (|2010l ') is 

i^LP = avCs/i/20 . (16) 

These relations are generated essentially by reversing the 
summation form of SPH equations which involve artificial 
viscosity back to continuum limits. Further details are given 
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in the respective papers. Note that the Gadget-2 definition 
of smoothing length has been used in each case. 



4 ESTIMATING LOCAL u FROM ENTROPY 
OR VISCOUS ENERGY PRODUCTION 

4.1 Motivation from fluid equations 

An artificial viscosity term (oc Hab) appears in both the mo- 
mentum and energy hydrodynamic equations. Each of these 
roles affects the evolution and structure of the system, as 
the trajectories of fluid elements are modified and as kinetic 
energy is transformed into heat. The standard form of the 
(isothermal) ring-spreading test has been used essentially to 
analyse the effective viscosity via the Euler equation, with 
the spreading being caused by the combined action of the 
viscous and several other additive term^, but it seems bet- 
ter to use, instead, the energy equation, where the effect of 
viscosity in generating entropy can be clearly singled out. 

Internal energy in fluid systems can be produced both 
by viscosity and by non-viscous compression. This is ex- 
pressed in the first law of thermodynamics, as applied to a 
co-moving fluid element with unit mass: 



de 



de = Tds - Pd{l/p) , 



(17) 



where T is the temperature and s is the specific entropy. 
The specific internal energy produced by viscosity, de^, cor- 
responds to the first term on the right hand side, and its 
time derivative is simply 



; Ts. 



(18) 



This quantity can be calculated directly from the SPH ex- 
pressions in all cases, whatever the flow geometry or equa- 
tion of state, giving a direct local measure of the effects of 
viscosity. The only caveat is that there should be no other 
sources of entropy production to contaminate the interpre- 
tation. 

For the particular case of the (widely used) class of 
polytropic equations of state implemented by Gadget-2, the 
derivative of Eq.|3] gives the following expression for de: 



de 



P 



,7-1 



-dK + Kp^-'^dp 



7-1 

£^dK^Pd(^- 
7-1 \p 



(19) 



where the relations P = Kp^ and d{l/p) = —dp/p^ have 
been used in arriving at the form of the second line. Match- 
ing terms with the flrst law in Eg. 1171 and referring to Eg. 1181 
leads to the following set of formulae linking K, s and e„: 



07-1 



(20) 



Within Gadget-2, computation of the evolution of ev is very 
similar to that for the entropic function, K (cf. Eqs. 1111 1181 
and l^ : 



^ In practice, ring-spreading is performed in 2D only, with the 
thermal pressure terms zeroed in a 'cold' disc approximation to 
isolate the viscous term, so that the 3D behaviour of the artificial 
viscosity is typically not measured. 



va 



7-1 



b 



(21) 



and this can be evaluated straightforwardly for each fluid el- 
ement. Note that in the further specialisation to an isother- 
mal equation of state with 7=1 (which we will be using 
below for making comparison with the ring-spreading test), 
then the entropic function, K, becomes a constant, but the 
true entropy, s, continues to change under the action of vis- 
cosity, as does e^. 

Importantly, e^ is determined from the contribution of 
only a single term which contains Uab- In quantifying the 
effect of any SPH artificial viscosity prescription on a disc 
system, this isolation greatly simplifies analysis (and inter- 
pretation). In general, entropy can be produced by shearing, 
normal compression and shock compression, the totality of 
which are accounted for when measuring ev, due to its gen- 
eral relation to specific entropy. In contrast, only the shear- 
ing is important under the assumptions of the ring-spreading 
test (or the analytic approximations). 



4.2 Local viscous heating of fluid elements 

Monitoring ev already provides a suitable way of demon- 
strating the local effects of artificial viscosity, and one could 
stop at that point. However, if one wants to calculate an ef- 
fective shear viscosity coefficient (for circumstances involv- 
ing shear and no significant compression), one needs to write 
an expression for ev in terms o f v and t he she aring velocity 
field (see, e.g.. Appendix B of iTassoul (|l978l )) and then to 
invert it to give an expression for f. For general planar rota- 
tional motion, one finds the following expression for specific 
internal energy creation ( Szuszkicwicz & Miller 1997): 



4 

ev = -. 



dVr 

dr 



dVr Vr 
J.2 Qj. J. 



+ V 



dr 



(22) 



using cylindrical polar coordinates (r, z) here and after. 
This form may be inverted easily for an expression of the 
effective viscosity: 



.(r) = 



3ev 



dVr 

dr 



+ 



dVr 

dr 



+ 



dr 



.(23) 



Each term on the right hand side is composed of variables 
which are directly obtainable from the SPH particles. As 
discussed above, the viscous heating is a direct measure of 
entropy creation due to viscosity, and the gradient terms 
may be estimated, e.g., within SPH particle kernels or from 
averaged ring values. 

In this study, we describe and utilise an expression 
which includes solely the term involving Q. This facilitates 
comparisons with existing approximations and, as shown be- 
low, leads to particularly simple expressions of u{r) (though, 
it must be emphasised, that the preceding expression is by 
no means prohibitively complicated). 

The surface area of a narrow, axisymmetric ring of 
width, Ar (centred around r), is 27rrAr, and the mass of 
the ring is therefore Am — 27rrEAr. Taking the matter to 
be moving on basically circular orbi ts, the viscous torque 
exerted on the ring can be written as l|Frank. King fc Raind 
I2OO2I) : 
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G(r) = — 27rrVE. 
dr 



(24) 



The work done by this torque leads to local heating in the 
rotating flow. The specific heating rate due to viscosity is 



dQ G{r) Ar _ f dn\ 2nr^uT,Ar 
dr Am \ dr / Am 



dO, 
dr 



(25) 



which is recognized as the final term in Eq. 1221 

After trivial rearrangement, the expression for the vis- 
cosity is given in terms of easily known quantities as: 

"(^)=("f ) (26) 

Compared to existing methods for evaluating effective vis- 
cosity in rotating flows, this expression is particularly gen- 
eral, having made no assumptions of rotation proflle, equa- 
tion of state, constancy of u or artiflcial viscosity prescrip- 
tion. Moreover, as a function of radius, this expression pro- 
vides a more detailed and useful description of viscous be- 
haviour throughout the simulation. As noted above, more 
general viscous stresses may be applied as well, such as that 
of Eq. 1221 In the following examples, we implement Eq. [^Hl 
with different disc stuctures and using different artiflcial vis- 
cosity prescriptions. 



5 EXAMPLES AND COMPARISONS 

5.1 Local heating applied to thin discs 

In the examples below, we apply the expression for the ef- 
fective local viscosity in Eq. [26] to test scenarios of ring- 
spreading and thin disc cases. Using assumptions of axisym- 
metry and circular, Keplerian orbits, the viscosity expression 
is 

Presuming vertical equilibrium with some thin disc scale 
height, H, the Shakura-Sunyaev a is then given as 



a{r) 



9 n|c,// 



(28) 



In the case of polytropic, nonisothermal equations of 
state, these expressions for viscosity and the Shakura- 
Sunyaev a may be written in terms of the entropic function 
as: 



v{r) 



9(7 



<r) = - 



9 (7 - l)^lc.H ' 



(29) 



(30) 



respectively. These last forms are particularly relevant to 
Gadget-2, where the entropic function is directly evolved, 
and they are used in obtaining the results presented below. 

Finally, it is interesting to note that, for the speciflc 
case of an isothermal EOS, Eq.[2S]may be expressed neatly 
as: 

«M = ^7^<x^, (31) 



using the expression for the scale height given in Eq. 1321 
below. Physically, in this latter form, a directly relates the 
local dynamical timescale (fdyn ~ ^~^) and the viscous ther- 
mal timescale (t th.v oc cj/f-v), as it should from accretion 



disc theory (e.g.. |Pringi3 (|l98ll ')') 



5.2 Thin Disc Structure 

We now briefly describe the structure equations for the two 
varieties of thin isothermal discs used in this study (in di- 
mensionless units). Each disc is taken to be in vertical, hy- 
drostatic equilibrium with respect to the gradient of the 
gravitational potential of the massive central object. The 
density distribution, scale height and central density (in the 
equatorial plane) are, respectively: 



p{r,z) = pc(r) exp 



H{r) 



H{r) = 



Pc(r) 



-I 1/2 



GMc 
E(r) 



Ok 



(32) 



[27r]i/2_ff(r-) ■ 

Since this Gaussian density proflle only asymptotically ap- 
proaches zero in the vertical direction, an arbitrary disc 
height for the boundary of SPH particles above/below the 
equatorial plane at any radius must be designated, and we 
set this to be three times the scale height in order to in- 
corporate > 99% of the mass of the theoretical distribution 
into the simulation (though, the non-arbitrary scale height, 
H , is used everywhere else in the analytic relations and dis- 
cussion) . 

The vertical velocity is taken to be uniformly zero, and 
the profile of the azimuthal velocity is very nearly Keple- 
rian, i.e., v^{r,z) « VK{r). The global sound speed, Cg, is a 
parameter restricted in accordance with the thinness condi- 



tion, H r 



< GMc/r, or w^/c^ 



1, where fi 



9f7K ci 



is essentially the Mach number of the (very supersonic) az- 
imuthal fiow. Using both the height from Eq.[32]and the def- 
inition of the viscosity parametrisation in Eq. 1131 the global 
disc viscosity can be related to the Shakura-Sunyaev a as a 
function of /i and the Keplerian specific angular momentum, 
Jk: 

rvK Jk .„„s 
u = a — — = a ^ . (33) 

Finally, the initial surface density of the flows must be 
supplied for use in Eqs. [321 The formulation of the spreading 
ring, Sr, has been given in Eg. 1141 We have also calculated a 
more disc-like system (of similar mass and surface density), 
with radial span at t = of (rmin, '"max) = (0.7, 1.4). Here, 
the arbitrary surface density proflle has been chosen to be 
flat in this region initially, T.F{r, 0) = Er(1, r = 0.01)/2, 
and the resulting disc mass was ~ 2. This is considered 
solely as a separate test case and not as a representation of 
a speciflc, physical system (e.g., the initial boundaries are 
too discontinuous for an equilibrium system). 

5.3 Global and local viscosity estimates 

Here, we flrst present results for ring-spreading tests with 
the Balsara correction included for various values of the 
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Table 1. Comparison of viscosity parameter values at multiple 
resolutions, from both local (entropic) and global (kinematic) 
methods for the same disc {fii = 30) shown in Fig.[l](at t = 2000). 



parameter 


method 


N = b X 10 


1 X 10^ 


2 X 10^ 


^avc 




0.058 


0.045 


0.036 


I'avo/ClO-S) 


entr. 


6.0 


4.6 


3.5 


I^E,m/(10-6) 


cntr. 


4.9 


3.7 


2.7 


v/(10-6) 


kin. 


3.5 


2.6 


2.1 




entr. 


0.056 


0.042 


0.032 




entr. 


0.046 


0.035 


0.026 


a 


kin. 


0.03 


0.02 


0.02 



1 

- Hi=30 


1 , , , , 1 , 


v=3.5xl0"' 


A : 


a=0.03 


|1 : 










^-.^gja^r,:.,^^ ^ 





parameter, /ii = fj,{r — 1), where r = 1 is the radius of 
the initial midpoint of the ring (and using units in which 
GMc = 0.01). First, the viscous properties are determined 
kinematically by using the global evolution of the surface 
density and fits to the theoretical progression given by the 
analytic formulation. Then, local values of v and a are cal- 
culated from iv using Eqs. I27II28I Unless stated otherwise, 
Ov = 0.8 in all of the presented simulations (and, in the 
notation of JJ^l the quadratic term coefficient, /3v ~ 2.4). 

The top panel of Fig. [T] shows the Er time evolution of 
a ring with /ii = 30 at a resolution of A'' = 5 x lO** (binned 
values plotted as diamonds). The value oi ly — 3.5 x 10~^ 
can then be inferred by matching with analytic fits given 
by Eq. 1141 which are also shown (solid lines) . Use of Eq. [33] 
then easily yields values of a « 0.03 at r ~ 1. The profile 
from the simulation roughly matches the analytic expecta- 
tion, and while fitting the curves, more weight was given 
to the locations of peak surface density, with the inner and 
outer margins being slightly over- and under-dense, respec- 
tively. It is striking that there is non-neglible effective shear 
viscosity present in these simulations, even with the Balsara 
correction. 

The local estimates of viscosity parameters coming from 
measuring K and evaluating Eqs. I27I28I are also shown 
in Fig. [T] where favc is a mass-weighted average between 
0.5 ^ r ^ 1.5, and !/E,m is the value for the annulus with 
maximum surface density (same notations for a). For a 
given flow, both sets of locally-estimated values of kine- 
matic viscosity and Shakura-Sunyaev a are roughly simi- 
lar to the globally-obtained curve-flts, though always having 
larger magnitudes. Comparing the results of the two meth- 
ods, the values at the peak surface density of the local cal- 
culations are much closer to those of the global estimates, 
which is not surprising given that the former emphasised 
the maximum-density regions in fltting. Results from higher 
resolution rings were similar but showed a slight resolution 
dependence in the effective viscosity values obtained (dis- 
cussed further below); a comparison of parameters is given 
in Table [U 

The profiles of the local v and a values (lower two pan- 
els of Fig. [1} are roughly constant with time (and also for 
the various resolutions tested), with minima near the loca- 
tion of maximal Er. This can be understood by examining 
the corresponding heating rate profiles in Fig. [2] In the ex- 
panding layers at r > 1, tv is nearly constant, and there- 
fore, by Eq. 1271 u oc r'^. The inner layers are compressed 
as they spread and have a much higher heating rate. In the 



+ v„, = 5.7x10-', Vj„ = 4.5x10', t = 1000. 
X v.„ = 6.0x10-', Vj„ = 4.9x10', t = 2000. 
o v„, = 6.2x10-', Vj„ = 5.2x10', t = 3000. 



+ = 0.052, a. 
K a^^e = 0.056, a. 



0.043, h,„ = 0.055 
0.046, h,„ = 0.058 
0.058, aj„ = 0.053, h.„ = 0.060 



Figure 1. (Here and following:) Top: standard u and a estimation 
from Er evolution ('o') with analytic solution (solid line), shown 
at t = 0, 400, 1000, 2000 and 3000. Middle: local u (average and 
value at maximum S) from Eg. 1271 Bottom: local a from v values 
in middle panel, using Eq. 1281 and average smoothing length, 
/lave. Af = 50k. 



regions containing the most mass (0.7 ^ r ^ 1-3), the values 
of viscosity and the Shakura-Sunyaev parameter are fairly 
constant, in accord with the overall model assumptions. 

The large viscosity values at the inner and outer disc 
boundaries (Fig. [l| correspond to very low density regions 
which are poorly resolved. A consequence of the smoothing 
length condition in Eq.[SJ used to stabilise the interpolations, 
is that particle kernels are artificially large at the edges of 
the simulation, and therefore the density is artificially low. 
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1.000 



0.100 



0.010 



0.001 




Figure 2. The (scaled) specific lieating, tv/Cg, at t = 1000 
2000 ('*') and 3000 ('o') from ttie same simulation shown in Fig.[T] 
The heating is mostly constant for the 'outflowing' material at 
r > 1.0 but is much higher for inner radii. 



Hence, Gadget-2 and other SPH formulations with similarly 
defined smoothing lengths necessarily require careful bound- 
ary conditions. As a function of time, the values of both Oavc 
and as.m, increase slightly (along with ftavc). 

Fig. [3] shows results from investigating other values 
of /ii (= 60 and 120), again using 5 x 10* particles. For 
both the global/kinematic and local/entropic measures, a 
increases with increasing /ii, whereas v decreases (with the 
local/entropic values being systematically higher than the 
global/kinematic ones in each case). The resolution trends 
are the same as those for ^\ = 30. It then appears that, 
while the value of the artificial viscosity parameter has re- 
mained constant (qv = 0.8 for all of these simulations), the 
kinematic viscosity and disc-a depend on both the sound 
speed and the resolution. To explain the former depen- 
dence, one can note that Cg appears directly in the expres- 
sion for Ilai, in Eq. |9l The latter dependence results from 
the fact that, inherently in SPH, the smoothing length of 
a fluid element determines the length scale of its interac- 
tions. With decreasing particle number, the cross-sectional 
spread of an element increases (/lave oc N^^^), resulting in 
greater shear in the flows. Using only t he bulk term i n Hat 
for studying non-axisymmetric discs, (|Murravl [l99^ ') esti- 
mated linear dependence for both sound speed and smooth- 
ing length (which was assumed to be constant), i.e., i/ oc Csh; 
recently similar analysis has been performed for the case 
of artificial tensor-vis cosity in studies of warped discs by 
iLodato fc Pricell2010l ). 

This dependency of effective viscosity on features which 
vary across the disc shows the importance of local measure- 
ment. Both the spatial and temporal variations in surface 
density away from the exact analytic solution have been 
noted above. In using the local estimates to determine pa- 
rameters with Eqs. 12711281 one has the dual advantage of 
avoiding any arbitrariness which arises in global curve- fitting 
and of characterising radial dependences quantitatively. This 
approach allows a more complete description of the be- 
haviour of the SPH artificial viscosity in the ring. 

For example, in the very thin /ii — 120 rings the effec- 



tive viscosity has a much smaller constant-value region than 
in rings with higher Cg. Profiles of a steeply increase to either 
side of r « 1, also making analytic curve-fitting less exact. 
Apart from differences due to compression and expansion of 
annuli, the curves reflect the influence of boundary condi- 
tions and the effectively poorer resolution within these flat- 
tened discs {H oc fJ.~^), even though the volume is smaller 
for the same A'^: since the ratio of surface area to volume is 
necessarily larger as well, the kernels of more SPH particles 
encompass empty space outside the disc, artificially increas- 
ing h (and similarly decreasing p). This would suggest the 
practical necessity for increased N in thinner, high-/i discs 
to avoid numerical effects in the viscous evolution and arti- 
ficially high accretion rates. In some sense, this is a state- 
ment of the intuitive notion that a resolution criterion such 
as H/h ^ (a few) is required for accurately characterizing 
discs. 

As a comparison, we now investigate rotating flows with 
flatter surface density proflles, for which the 'empty kernel' 
boundary behaviour should have less effect on the evolving 
system. We therefore use local viscosity measures to analyse 
the systems with the Ep deflned in ij5.2l which are initially 
more 'disc-like' in the sense of having a greater radial extent 
of mass. Time-evolution proflles of surface density are shown 
in Fig. Ufor /ii — 60 and /ii = 120. Indeed, for these discs 
the values of i/ and a are fairly constant functions of radius, 
and, for a given pi, the measured viscosity values agree well 
with those from the ring tests. 

In order to quantify the behaviour of the Balsara- 
corrected artiflcial viscosity prescription, 11^^ with a given 
Qv, for rotating flows with different resolutions and sound 
speeds, the resulting Shakura-Sunyaev parameter values for 
all rings and discs from this study were plotted in /lavo-Qave 
and ftavc-Q!E,m planes (Fig. [S]). Note that the a values ap- 
pear to vary quite linearly with smoothing length, with pi 
determining the gradient; the values of as.m appear to be in- 
dependent of the specific profile of surface density, whereas 
disc values have noticeably lower values of Oave than the 
rings. This (small) qualitative difference is to be expected, 
as boundaries affect the former less. While the values in both 
plots appear to be quasi-convergent towards a = as ftave 
tends to zero, it is not possible to be precise about this due 
to the approximate nature of the plot. 



5.4 Further examples and altering artificial 
viscosity 

5.4-1 Evaluating the artificial viscosity in 
Balsara- corrected flows 

The simulations presented above have all utilised the same 
value of Qfv = 0.8 and the Balsara-corrected form for the 
artificial viscosity, Il'ab- In this section we investigate the ef- 
fects of varying Qv, using only the local entropy production 
method to measure the viscous properties of the fiows. It 
should be noted that for these predominantly smooth and 
non-shocking flows, the linear term of the artificial viscos- 
ity dominates entropy production. Neglecting the quadratic 
term (/3v = 0) produces very little change in the disc fiow. 
Therefore, we only further consider here the effects of vary- 
ing Qv and of using or not using the Balsara correction. 
With n^i,, even significantly different values of Qv result 
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Figure 3. Ring-spreading test with A'^=50k, shown at times given in the middle panel (and at i = 0, top). Left, fii = 60; right, fii = 120 
(cf. Fig. [T] for notations). 



in only minor changes of the disc parameters, as was to be 
expected. Originally, the ring with fii = 30 and — 0.8 
(seen in Fig. [T]) had (aavc, QE.m) values of approximately 
(0.056, 0.046); increasing to large values of = 1.6 and 2.4 
resulted in (0.060, 0.049) and (0.066, 0.057), respectively; 
only at extremely high values were significant changes seen, 
e.g. av = 8.0 yielded (0.109, 0.091). Thus, in these purely 
rotational flows, the shear-correction does limit the artifi- 
cial viscosity term (although we stress again that it remains 
significantly non-zero), and the values shown in Fig. [5] are 
roughly independent of within the range studied. From 
fitting various disc parameters for results from the preced- 



ing simulations (across a range of radii, 0.5 < r < 1.5), it 
emerges that the effective viscosity is well represented by 
the following (simple) expression: 



a{r) ^ - 



hir) 



2H{t) 



(34) 



explicitly noting that a generally varies with position, as 
both the smoothing length and disc height are functions of 
radius. Therefore, when the Balsara correction is applied, 
the interaction length among particles appears to strongly 
determine the viscous dissipation in the flow. 
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= 0.083, ay „ = 0.063, h,„ = 0.046 

= 0.086, ay „ = 0.072, h.„ = 0.047 

= 0.087, ay „ = 0.066, h,„ = 0.047 

= 0.089, ar„ = 0.080, h.„ = 0.048 



= 0.133, aj„ 
= 0.156, ay „ 
= 0.154, ay „ 
= 0.156, ar„ 



= 0.175, h,„ = 0.036 
= 0.097, h,„ = 0.037 
= 0.107, h,„ = 0.037 
= 0.112, h.„ = 0.038 



Figure 4. Flat-profile discs with constant 12p{x,t = 0) (as described in i|5.2ll and N = 5 X 10^. Left, fii = 60; right, /ii = 120. 



5.4-2 Evaluating the uncorrected artifical viscosity inflows 



In switching from Il'^i, to Hat (no Balsara correction), the 
eflfective viscosity for Qv — 0.8 is nearly doubled. The clearly 
linear dependence of aE.m and aave on in this case (and 
with /3v = Sqv) is shown for a series of discs in Fig. [G] Best 
fit lines of the data points are provided as visual guides, with 
open and skeletal symbols denoting Ofavo and a^.m, respec- 
tively. The dependences on both Cs and ftavo (which is mainly 
determined by A'^) remain qualitatively similar to those ob- 
served in the Balsara corrected simulations (i.e., increasing 
either parameter leads to an increase in q). In the range of 
these examples, the 'empirical' dependence of the Shakura- 
Sunyaev a on SPH and How parameters can be derived from 
a fitting procedure similar to that used for Ea. l34l above and 
given as: 



Q(r) ^ 0A3a^fi{ry^^^ 



h{r) 



-1 3/2 



2H{r) 



(35) 



The functional form offers direct insight into the depen- 
dences of the effective viscosity on disc parameters, and this 
relation can be compared with existing approximate, ana- 
lytic formulae in the literature. Fig. [7] shows a co mparison 
of V d erive d from Eq. | 35 | w i th for mulae derived by iMurravi 
(l99d) and iLodato fc Pried (|2010|). Even though the latter 
two relations were applied for different SPH schemes (as well 
as containing some dependence on both differing kernels and 
dimensionality), the figure shows that there are regions of 
good agreement with the empirical values, namely for small 
smoothing lengths. (We note that the curves derived in the 
present study have been restricted to the range of h val- 
ues observed in simulations.) However, as h increases, the 
differences from the analytic relations increase significantly. 
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Figure 5. Phase space of ciE,m (left) and Oavc (right) vs /lavc from all of the rotating flow simulations in this study using n^j^ and 
Ov = 0.8. Ring and flat-disc simulations (Sj^ and Sp) are represented, respectively, by 'o' and '+' for fii = 30; by '□' and 'x' for 
fil = 60; and by 'A' and '*' for fii = 120. The points appear aligned by azimuthal mach number, and best fit lines for /ii = 30, 60, 120 
are included as visual guides. 
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Figure 6. Phase space of Shakura-Sunyaev alpha values vs SPH 
artificial viscosity, using the non-Balsara corrected W-ab- For Ov = 



0.2, 0.8, 1.6, 



values are respectively shown in open symbols 



(A, □ and o) and aj; ^ values, in skeletal symbols ('*', '+' and 
' X '). The different values of ii\ and /lave (the latter mainly varying 
with A'^) are shown on the figure, with best fit lines of the data 
points as visual guides. 



Formulae such as Eq. [35] can be useful in efficiently es- 
tablishing the setup of a simulation (e.g., choosing the num- 
ber of particles to produce a required smoothing length, as 
a sort of resolution requirement). For the non-Balsara flows, 
extrapolation far beyond the examined parameter space 




Figure 7. A comparison of analj'tic expressions for t he effec- 
tive d isc viscosit y, v. Previous an a lytic expressions from [Murray! 
lll996l ) and from iLodato fc Pricel 1 I2OIOI) are shown with dotted 
and dashed lines, respectively; solid lines show Eq. 1351 for repre- 
sentative radial values. In this ring simulation with TV = 5 X 10*, 
the h values were mainly between 3 X 10~^ and 2 X 10~^, hence 
the displayed interval of the Eq. 1351 curves. 



with any approximated function is most likely unsuitable. 
However, a function such as this coming from tests with a 
fairly wide range of values for the independent quantities 
may provide an adequate level of accuracy for interpolation. 



6 SUMMARY AND DISCUSSION 

An entropy-based local method has been described to mea- 
sure the effective kinematic viscosity in the flow due to 
SPH artificial viscosity, which can be simply related to 
the Shakura-Sunyaev a for appropriate thin discs. This ap- 
proach has several advantages over that of estimating a 
global V from an isothermal ring-spreading test. First, it 



© 2011 RAS, MNRAS 000, [T]-?? 



12 P. A. Taylor and J. C. Miller 



calculates the local viscosity of an annulus from the (vis- 
cous) energy or entropy expression, which contains only the 
term dependent on Ilab, and therefore only viscous effects 
are measured. This procedure then yields f and a as func- 
tions of radius, so that variations in the effective values 
across a disc can be determined; global values may be ob- 
tained easily via averaging, if desired. In addition, results 
are generated immediately, without the necessity of finding 
fits for analytic curves as in ring-spreading, and also thereby 
eliminating arbitrariness from the estimates. Moreover, the 
viscosity values are calculated from quantities which are al- 
ready calculated in the SPH algorithm. 

Finally, whereas traditional ring-spreading tests are re- 
stricted to use with an isothermal equation of state (or, in 
practice, even more restricted to cold, 2D systems) because 
it is only then that one has the analytic solution for com- 
parison, this local method can be used with any equation of 
state. The effective viscosity may be obtained for rotating 
flows in general as a function of viscous dissipation, such 
as from Eq. I26b r 1231 (both derived using the very general 
Eq. I18|l . The exact form would simply differ from that of 
Eq. 1271 in which a Keplerian profile (SI — SIk) has been 
assumed. Further work will also investigate an analogous 
approach applicable to estimating the effective local shear 
viscosity for systems in linear motion. 

We have compared two existing methods (standard 
ring-spreading and analytic approximations) with the new, 
local entropy method in studying the behaviour of the SPIf 
artificial viscosity in rotating fiows. In cases where the par- 
ticular assumptions of the former methods were applicable, 
average results of all approaches were shown generally to 
be in fair agreement, with far greater detail coming from 
the local entropy method. Further examples changing the 
disc structure and the SPH artificial viscosity values were 
given. Examples of explicit expressions which approximate 
the dependence of the numerical a on various fiow and SPH 
parameters were given for both the Balsara and non-Balsara 
corrected forms of artificial viscosity. Procedures such as 
these may provide useful guides when setting up simulations. 
However, further work must be done to test the generaliz- 
ability of these specific forms, i.e. Eqs. I34I35I across wider 
parameter space and to other SPH codes. The means for 
estimating u and a in each case, however, applies generally 
with the local method. 

We note that the relations between SPH variables and 
the effective viscosity in disc flows described here can serve 
as useful examples of the type of analysis possible with 
the local entropy measures. The ability for hydrodynamic 
schemes to represent flows with particular physical viscosity 
characterizations, such as the widely-used Shakura-Sunyaev 
discs, may be investigated more directly and in more detail 
with the local method. Moreover, this analysis may be use- 
ful in the further development and evaluation of physically- 
motivated viscosity terms in SPH, to which existing meth- 
ods cannot typically be applied. For example, it was shown 
here that the a in the case of Balsara corrected artificial vis- 
cosity remained finite and linearly dependent on smoothing 
length. Studies continue to develop schemes for minimizing 
unwanted effects of ar tificial viscosity in rotating flows, such 
as the recent work bv lCullen fc Dehneni ({201Q '). and we sug- 
gest that the simple approach outlined here for measuring 
the local viscous effects directly would be extremely useful 



in quantitatively assessing these developments. (Indeed, we 
will return to discussion of this in a following paper.) 

Finally, several studies have described and in- 
vestigated numerical implementations of physical vis- 
cosity prescriptions derived direc tly from th e ana - 
lytic Navier-Stokes equations (e .g., Flebbe et al.l l|l994 ): 
iLanzafame. Belvedere fc Moltenil (|2006l )). with separate 
bulk and shear viscosity terms (unlike the standard SPH ar- 
tificial viscosity, which does not differentiate the two); this 
is what one must do in order to incorporate real viscos- 
ity behaviour into SPH calculations. However, some numer- 
ical considerations remain for such physical viscosity pre- 
scriptions, in particular regarding spatia l derivatives and de- 
pendence on particle disorder (see, e.g.. IRosswo 3 (|2009l). as 
well as Appendix 1X1 of the present paper). While one could 
think of using the SPH artificial viscosity as a model for a 
real, physical one (see , e.g., the discussion in Sec. 3.2.3. of 
iLodato fc Pried (|2010l )). this provides only a rough approx- 
imation, and it is best to think of the present work in terms 
of establishing acceptable limits for the setup parameters so 
that the artificial viscosity does not interfere with the effects 
of the real one. 
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APPENDIX A: INITIAL PARTICLE 
DISTRIBUTION 

As any SPH simulation evolves, the particles asymptotically 
tend toward an equilibrium of mutual separation distances 

i i.e., no artificial clumps or voids) within the local profile 
MonaghanI 120061 ). The presence of either lattice-regularity 
or purely random (over)density fluctuations in the initial 
conditions impedes this process, resulting instead in purely 
numerical features, the effects of whi ch may continue to 
propa g ate in the simula ted flows, e.g., Ilmaeda fc Inutsukal 
l|2002h : lMonaghanl (| 20061 ). For example, locating SPH parti- 
cles at regular gridpoints often leads to spurious spiral modes 
in rotating systems (though 'real' spirals may form due to 
su sceptibili t y con d itions in certain cases, such as those given 
bv lToomrd l^l964^ : ISpeith fc Kiev! (|2003l )). In particle setup 
methods using random placements of equal-mass particles, 
some of these numerical side effects are limited, but fluctu- 
ations may still occur in the flow which are of an entirely 
artiflcial nature, particularly due to localised overdensities; 
typical solutions of these difflculties have been to enforce a 
minimum initial separation between randomly placed par- 
ticles (a computationally expensive, A'^'^ procedure which 
is not guaranteed to remove all fluctuations), as well as 
spending time at the beginning of a simulation to evolve 
the system in the presence of strong damping, while main- 
taining the desired profiles of density, velocity, energy, etc. 



jRosswog. Speith fc Wvnnll2004l : iMonaghanI [20061 : |r osswogi 
l2009h . 

In this work, we have utilised a new method to create 
equal-mass particle distributions which are both smooth and 
random (SAR) from the start of a simulation, i.e., limiting 
numerical fluctuations due to initial placement. The SAR 
method produces such conditions a.t t — 0, without the need 
for special runtime considerations (such as damping, etc.) 
during the simulation. There is very little fixed cost in cre- 
ating such a distribution, and neither damping nor further 
non-hydrodynamic numerics are introduced into the model. 
The SAR method for setting up initial conditions is valid 
for any number of dimensions and does not assume any par- 
ticular symmetries. 



Al SAR process steps 

To make SAR initial conditions for a simulation, an inter- 
mediate distribution of point particles is first created, called 
a 'glass'. The glass is defined by the property that its con- 
stituent particles are equidistant from their nearest neigh- 
bours but without regular structure (such as in a crystal). 
For practical considerations of generating a glass, it is com- 
putationally efficient to first make a single, small glass dis- 
tribution and then to tile copies with periodic boundaries 
appropriately matched to preserve interparticle spacing. 

The general process of creating an SAR particle distri- 
bution may be divided into the following steps: 

(i) seed a mini-glass distribution: place a small number, 
A'g, of point particles randomly in a cube of normalised edge 
length, Lg = 1, with periodic boundaries 

(ii) create the mini-glass: evolve the particles with mu- 
tually repulsive (here, inverse-square) forces, asymptotically 
approaching a constant density of equidistant particles (par- 
ticle velocity is zeroed after each timestep to limit 'over- 
shooting' equilibrium) 

(iii) create full glass: tile copies of the mini-glass cubes 
(into any given size/shape), with the periodic boundaries of 
tiles matched to preserve particle spacing 

(iv) fill simulation model: select a small volume of the 
glass (including particles) and map it into the model, ad- 
justing the element volume (and relative particle locations 
therein) to the local density of the profile by compressing or 
stretching it 

(v) repeat the previous step, smoothly mapping neigh- 
bouring volumes from the glass to neighbouring volumes in 
the model profile, maintaining an absence of spurious parti- 
cle number density fluctuations. 

Using small volumes, such a process preserves the relative 
spacing of all points while creating an appropriate density 
profile. 



A2 Testing the SAR distributions 

Here, we first tested a SAR distribution model using a 
standard SPH test of particle resolution and convergence: 
the collapse of a self-gravitating, isothermal sphere ini- 
tially in solid body rotation and with a non-axisymmetric 
(m — 2) density perturbation. The system fragments, 
and this has been well-examined using both Lagrangian 
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Figure Al. I sothermal sphere test ( equatorial) density contours with four density contours are drawn per decade. Col. 1: 'solution' 
from Fig. 4 of lBate fc Burkerj ^99^, for which = 8 X 10"'; peak densities at t = 1.24 and 1.26 are pmax/po ~ 2,400 and 78,500, 
respectively. Cols. 2 and 3: NEM and SAR distribution simulations, respectively; in both, N = 6.7 X 10'*. SAR results show much greater 
agreement with the first column 'solution', both qualitatively and quantitatively. 



and various Eulerian grid codes ("Boss fc BodcnhcimeJ I 19791 : 
iBate fc BurkertI 1 19971 : IXruelove et al.. 1998). SAR models 
were tested with total particle number A'' — 3.4 x 10*, 
6.7 X 10* and 2x 10^ and results were compared with lattice- 
like, nonequal mas s (NEM) place ments using a standard 
configuration (e.g. 'SpringeJ (|2005h ) at similar resolutions. 
Results of both simulation s were also compared w ith con- 
vergent 'solution' results of iBate fc BurkertI (|l997l ). 

Briefly, the initial conditions of the test are: at t = 0, 
we start with an isothermal [P — c^p) sphere having uni- 
form sound speed, Cg — 1.66 x 10* cm s~^; radius, Rs — 
5 X 10^® cm; mass, M — IM©; and solid body angular ve- 
locity, Lj — 7.2 X 10~*'^ rad s~^. The density distribution has 
an m = 2 perturbation over the azimuthal angle given by 
p{<j}) = po[l + 0.1 cos(2(^)], where po = 3.82 x 10"^* g cm'^ 
All quantities in this section are scaled by the following char- 
acteristic quantities: the initial value of the outer radius, Rs; 
the free-fall time, tff = {3Tv/32Gpo)^^^; and the mass, Mq. 

Fig. lAll shows plots of density in the equatorial plane 
of the SAR and NEM distributions with = 6.7 x 10* par- 
ticles, as well as the convergent 'solution' (top), a,t t = 1.24 
and 1.26. Quantitatively, the peak density of the SAR sim- 
ulation was much closer to the expected value than that 
of the NEM case (and at fairly low resolution). This re- 
mained true for all N values investigated, and furthermore, 
the SAR contours were significantly smoother and more 'flu- 
idlike' with fewer numerical fluctuations. The simulation re- 
sults using SAR conditions converged much more quickly to 
the expected results as resolution was increased, an impor- 
tant property for creating efficient simulations. 



As an additional example of SAR results. Fig. IA2l shows 
images of particle locations for an evolving SAR ring sim- 
ulation (the same system as in the left panel of Fig. [3]), 
projected onto the equatorial plane. Over time, the particle 
distribution within system retains its good qualitative prop- 
erties of smoothness and randomness. Note the absence of 
any temporary artificial spirals, clumps or sub-structures, 
such as often appear in SPH simulations, even at very early 
times. Both the efficiency of modelling and the relative ease 
of formulation (the computation requires only a few simple 
loops through N particles), in addition to a reduction of 
numerical artefacts, make the SAR initial conditions very 
beneficial for use in SPH simulations. 

Importantly, the SAR particle method is effective be- 
cause the quasi-continuous mapping retains the relative 
spacing between particles from the initial glas^. The new 
density profile is obtained smoothly and accurately, with- 
out the presence of fluctuations or numerical artefacts. How 
smoothly a density profile is reproduced is principally de- 
termined by the size of the incremental volume element. 
The SAR method does result in a small amount of localised 
asymmetry in the particle distribution, as volume elements 
are not generally scaled uniformly in all directions, but for 

^ A version of particle set up obtained by warping a constant den- 
sity profile llHeranj|l993 'l has been pointed out by a referee. Al- 
though the mapping feature therein is ostensibly similar to what 
is described here, that method was limited for use with a cer- 
tain subset of spherically symmetric profiles and was made from 
uniformly placed (not SAR) particles. 
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t=40 




t=1200 




t=7600 





Figure A2. A sequence of positions of SPH particles, projected onto the equatorial plane, for the fii = 60, TV = 5 X 10^ model shown 
in the left column of Fig. |3] Note the early and continued absence of clumps, voids and temporary spiral. 



small elements this did not appear to produce spurious arte- 
facts during simulations. Improvements to the method are 
currently being considered. 
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